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Abstract 

Performing electronic structure calculations for large systems, such 
as nanoparticles or metal clusters, via orbital based Hartree-Fock or 
Kohn-Sham theories is computationally demanding. To study such 
systems, therefore, we have taken recourse to the hydrodynamic ap- 
proach to time-dependent density-functional theory. In this paper we 
develop variation-perturbation method within this theory in terms of 
the particle and current densities of a system. We then apply this 
to study the linear and nonlinear response properties of alkali metal 
clusters within the spherical jellium background model. 
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I Introduction 



The hydro dynamic analogy of quantum mechanics was first explored by 
Madelung [l[] who transformed the single particle Schrodinger equation into 
a pair of hydrodynamical equations. The theory views |2], || the electron 
cloud as a classical fluid moving under the action of classical Coulomb forces 
augmented by the forces of quantum origin. The basic dynamical variables of 
this theory are the particle and current densities which satisfy two fluid dy- 
namical equations, namely, the continuity and an Euler type equation. The 
work of Madelung was followed by Bloch's attempt [|J to develop hydrody- 
namical theory for many-electron systems within the realm of Thomas- Fermi 
(TF) theory 0, [/J. Although then proposed without any rigorous founda- 
tion, the theory can now be derived || from the equations of time-dependent 
density-functional theory (TDDFT) ||, 10|. It is based on the assumption 
that the dynamics of many-electron system can be described by considering 
it as a fluid of density p{r,t) and a velocity field v(r, t) which is assumed 
to be curl free (that is v(r,t) = — V5*(r, t), where S(r,t) is scalar velocity 
potential). 

Using p(r, t) and S(r, t) as conjugate variables Bloch derived two fluid 
dynamical equations: the continuity equation 



^ + V-(pVS) = 0, (1) 



and the Euler equation 



dS 5T . . r p(r't) , , 

Here T is the TF kinetic energy (KE) functional and v ext (r, t) represents the 
external potential. These equations were subsequently used to study pho- 
toabsorption cross section and collective excitation of atoms|jll|, collective 
excitations [12] and plasmons of metal clusters |13| and surface plasmons 



14, 15 in metals. 
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Bloch's theory also formed the basis of initial attempts by Ying [ITJ to 
extend density-functional theory (DFT) to include time-dependent (TD) ex- 
ternal potentials. He did this by replacing the TF KE term To by a general 
functional G[p] consisting of the KE and the exchange-correlation (XC) con- 
tribution to the total energy. In Ying's work it is implicit that, like in the 
static DFT, a universal functional G[p(r, t)] can be written for the TD prob- 
lem. The ad-hoc nature of Bloch's theory and its extension by Ying was 
removed by the pioneering works of Deb and Ghosh ||17|| , Bartolotti [18 and 
Runge and Gross Runge and Gross rigorously proved the existence of 
a Hohenberg-Kohn |20j like theorem for TD potentials, and showed that the 
TD density p(r, t) can be determined by solving the hydrodynamical equa- 
tions 

| + V.j = 0, (3) 
which is the continuity equation, and the Euler's equation 

§ = P[>(r J t)]. (4) 

Here P is the three-component density-functional of Runge and Gross and 
the vector j is the current density corresponding to the many-body wave- 
function \I/(ri, ■ • • , rjv; £). An explicit expression for P[p(r,t)] in terms of 
the wavefunction has recently been given pi using the TD differential virial 
theorem [Ell . 



Although there have been some calculations, as mentioned above, in the 
past by employing hydrodynamic theory, its full potential remains unex- 
plored. This is evidently because with the increasing computing resources 
one can perform |22j orbital based calculations like TD Hartee-Fock (HF) or 
TD Kohn-Sham (KS) with relative ease. Recently, however, hydrodynamical 
theory is being applied in situations where such orbital based calculations 
are still computationally difficult to implement. One such example is sys- 
tems which contain thousands of atom such as nanoparticles and clusters. In 
these systems hydrodynamical theory becomes the method of choice. Thus 
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the theory has been applied to study photoabsorption cross-section of metal 
particles P3fl , collective |23| and magnetoplasmon excitations of confined 
electronic systems and also to study the interaction of strong laser light with 
atomic systems PB|, ^7], ffiifl . Besides the computational ease offered by it, 
one is also tempted to work within the hydrodynamic formulation because it 
provides an intuitively appealing approach to the time dependent many-body 
problem. 

In this paper we develop perturbation theory within the hydrodynamic 
formalism to calculate linear and nonlinear response properties of large sys- 
tems. The motivation for this comes from our experience with the calculation 
of static response properties |29|, |3(J employing density based perturbation 
theory [[31] within the Hohenberg-Kohn formalism of DFT. The density based 
method reduces the numerical effort required for such calculations substan- 
tially while leading to reasonably accurate results |29], [3(J for the response 
properties. In the same manner, the hydrodynamical approach proves to be 
useful for calculating frequency dependent response properties of extended 
systems for which orbital based theories become quite difficult to implement 
because of the large number of orbitals involved. 

The work presented in this paper is divided into two parts: First, we de- 
rive the generalized Bloch type equation using the concept of time-averaged 
energy (quasi-energy) of an electronic system subjected to a TD periodic 
field. For calculating optical response properties we then develop variation- 
perturbation (VP) theory in terms of the quasi-energy using particle and 
current densities as the basic variables. This is presented in section II. The 
perturbation theory developed here proceeds along the lines of density based 
stationary-state perturbation theory |3l| by making use of the stationary 
nature of the time-averaged energy with respect to p(r,t) and S(r,t). In 
the second part we demonstrate the applicability of the perturbation the- 
ory developed here by calculating frequency dependent linear and nonlinear 
polarizabilities of inert gas atoms (section III) and comparing our numbers 
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with the standard results obtained from the wavefunctional approach. Hav- 
ing demonstrated the accuracy of hydrodynamical approach, we then apply 
it to calculate frequency dependent response properties of metal clusters with 
number of atoms up to 5000 using the spherical jellium background model 
(SJBM). 

II Variation perturbation method in hydrody- 
namical theory 

A. Time-averaged energy 

The central quantity around which the VP theory is developed in the static 
case is the ground-state energy of the system. For periodic TD hamiltonians, 
this role is played by the time-averaged energy or quasi-energy B3]. There- 
fore, in the following we first derive an expression for the quasi-energy as a 
functional of particle and current densities and show that it obeys stationar- 
ity for the correct solutions of these functions. As expected from the work of 
Runge and Gross [[0|], stationarity with respect to the density leads to the 



equation of motion for the density. In addition, variation with respect to the 
current density gives the continuity equation. 

We begin with the TD Schrodinger equation for the many-body wave 
function \l/(ri, • • • , r^; t) given by 

(tf(t)-^)*(r,t) = (5) 

where 

H(t) = f + V ee + V(t). (6) 

In the above equation T and V ee are the kinetic and electron-electron in- 
teraction energy operators, respectively, and V(t) denotes the TD external 
potential containing both the nuclear and the applied potential. For periodic 
hamiltonians, that is 

H(t + T) = H(t), (7) 
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where T is the time period, in accordance with Floquet's theorem there exists 
a solution \l/(r, t) of the form 

tf(r,t) = §(r,t)e- tEt (8) 

where $(r x , • • • , r^; i) is also periodic in time with time period T, i.e., $(£ + 
T) = $(t). Such a state has been termed as the steady-state of the system 
with E being the corresponding quasi-energy. The equation of motion for $ 
is easily seen to be 

\H(t)-i^$(r,t) = E$(r,t). (9) 

The corresponding expression for the quasi-energy is the time averaged ex- 
pectation value 

El<S>] = l [ (<S>\H(t)- i ^ t \<S>) ] j. (10) 

The curly bracket in Eq.fllOD denotes the time averaging over one period T 
defined as 

{fg} = f < [nt)9{t)dt (ii) 

The quasi-energy represents the average energy of induction |33j of a system 
subjected to a TD potential as is easily seen by the TD Hellmann-Feynman 
theorem [0 . 

To convert Eq . (|10D into its hydrodynamical counterpart we decompose 
the complex steady-state many-body wavefunction in polar form, so that 

$M) = X M)e^ (12) 

where both \ and S are real functions of r 1; r 2 , • • ■ , and are periodic in 
time with time period T. Further, S(r, t) also has a purely TD component 
which integrates to zero over time period T (for detail see Ref. |33[]). Note 
that, S is zero for the ground-state of the system. By substituting Eq.([L2|) 
in Eq. (|T0|) , the expression for the average energy becomes 

E[ X , S] = { ( X \f' + V ee + f^lx) - ^(X\^) } (13) 
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where 



^ = E(-^ 2 + ^|v ? ,5| 2 ) 



(14) 



Since the periodicity and reality of x(r, t) implies that 



= o, 



(15) 



the quasi-energy is given as 




' 2 + V 



ee 




Now by invoking the Runge-Gross theorem ||19|| , it can be written as a func- 
tional of the density alone. However, in the hydrodynamical formulation 
the density and the velocity potential S are treated as independent variables 
which means that the energy above is a functional of these two quantities. 
An advantage of this decoupling of the density and S is that one does not 
have to know the functional dependence of the S on the density. Further, this 
facilitates approximating the expectation value (x\ — |V 2 |x) as a functional 
of the density by the KE functionals well studied in static DFT. Evidently, 
the first three terms of the equation above can be represented by a functional 
of TD density as 



where fo( r ) represents the static external potential and TD part of the po- 
tential is represented by v app (r, t). This is because changes in S do not affect 
their values. The universal functional {F[p(r, t)]} given by 



{F[p(r,t)]} = {T s [p(r,t))} + {E H [p(r,t))} + {E x Mr,t))}- (18) 



Here {T s [p(r, t)]}, {^^^(r, t)}} and {^^^(r, t)}} represent the time-averaged 
KE, Hartree energy and exchange and correlation (XC) energy functionals 




(17) 
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respectively for the TD system. The TD particle density is given by 

P(r, t) = J x*(r, r a , • • • , r N , t)x(r, r 2 , • • • , r^, t)dr 2 • • • drjv (19) 

So far we have written the first three terms in terms of the particle density, a 
quantity defined in 3D configuration space. The last two terms representing 
the current still have all the co-ordinates of the configuration space in them. 
As such any equation involving 5(r 1; ■ • • , r^; t) can not be projected on to 
3D space. To do this one needs to consider some approximate form for the 
phase S. One such approximation for S which is generally employed [|3| is 
that it can be written as the sum of single particle phases, that is 

N 

S(r 1 ,---,v N ;t)=Y,S(r t ;t) (20) 
i=i 

with the same function S representing each electron. This approximation is 
equivalent to assuming the velocity field of the electron fluid to be curl free 
as was done by Bloch || in deriving Eq.(Q). With this approximation the 
average energy functional of Eq. (|i~6[) is given as 

E[p,S] = ^F\p(r,t)}+J{v (r)+v app (r,t))p(r,t)dr 

+ i|p(r,0(V5) 2 dr + /^p(r,0rfr| (21) 

This is the expression for the quasi-energy of a many electron system (under 
the approximation made above) interacting with a TD periodic potential. 
Since the purely TD component of S(r, t) is the same as the phase of the 
wavefunction, it does not contribute to the energy. In Eq.|2T| this is ensured 
by p(r, t) integrating to a fixed number of electrons at all times. We therefore 
drop it and work with only the co-ordinate dependent component of S(r,t). 
This is similar to separating out the overall phase of TD wavefunction |33|, |32| 
in the TD perturbation theory. 

We now demonstrate the variational nature of E[p, S] with respect to p 
and S. Making E[p, s] stationary with respect to p and S gives the Euler 
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equation 

K*) = + \(VS) 2 + v (r) + v app (r, t) + ^ (22) 

where p(t) is the Lagrange- multiplier ensuring that p(r, t) integrates to the 
correct number of electrons at each instant of time, and the continuity equa- 
tion 

^ + V ■ (pVS) = 0, (23) 



respectively. Eq.fl22|) is the same as that proposed by Ying|T6[. As such if 
F[p] is approximated by the TF functional, it gives the Bloch's hydrody- 
namical equation correctly. Further, for time independent hamiltonians it 
correctly reduces to the Euler equation of static DFT. All these facts demon- 
strate the variational nature of E[p, S] with respect to p and S. Employing 
this we now develop the VP method in terms of the particle and the current 
densities. We show that the (2n + 1) theorem and its variational corollary is 
satisfied in terms of these variables. 

B. Perturbation theory 

To develop perturbation theory we assume that v ex t(r, t) is relatively weak 
and under its action the ground-state density p^°'(v) changes to p(°'{r) + 
Ap(r, t) and the velocity potential changes to S^°'(r)+AS(r, t). The particle 
density change Ap satisfies the normalization condition 

[ Ap(r,t)dr = 0. (24) 

However no such condition is required for the change in the velocity potential 
AS. The changes Ap and AS* are expanded in perturbation series as 

Ap = £>Cfl 

3 

AS = J2 si3) i ( 25 ) 



9 



where p^' and correspond to the jth order terms in the perturbation 
parameter. The energy corresponding to + Ap and + AS is given by 



E\p<® + Ap, S<® + AS] = |F[p(°) + Ap]+J {v + v app ){p^ + Ap)dr 

+ - J V(S (0) + AS) ■ V(S (0) + AS)(p (0) + Ap)rfr 
+ /^±^»> + Ap )<ir } ,26) 

Using Eq.(EBp we now obtain the energy changes to different orders in per- 
turbation parameter employing an approach identical to the one adopted in 
Ref. |U| for time independent density based perturbation theory. The result- 



ing expressions for average energies to different orders are: 

E {1) = {/««(r,ty°>(r)dr}, (27) 



Em - {y wS,.) A, ' tw '' ) ^ + /' ar '^ i(# 

+ / gg^M pW (r, t)dr + i |(VS« • VS« )p<°» (r)drl , (28) 
+ |(VS (1) - VS (1) )p (1) (r,t)c/r|, (29) 

+ a / mM&, m*. r/ 1 ' (r - (r ' t)pl " (r "- ' ) '" ) (r 



x drdr'dr"dr'" 
+ 



- J (VS™ . V5( 2 ))p(°)(r)rfr + - J (VS« ■ VS^)p^(r,t)di 
+ J(VSU-VSW)pW(r,t)dr + J ^P (2) (r, t)drj 



(30) 



and 

1 r 5 3 F[p] 



E{s> = {5/ ^,, t )^,V.^ )(r - f) ^ ( ^^''^- f)<M,,<y 

+ S / MM)M..wi-..)M^^ '" (r ' i)pU,(r ' i)pa,(r '' t) " W(r "' 



x drdr'dr"dr"' 

+ — f o w (r t)o {1) (r' t)o {1) (r" t) 

+ 120 7 5p(r,t)5p(r',t)5p(r",t))Sp(r'",t)5p(r"",t) P { , )P 1 , )P 1 , J 

x p« (r w , t)p^ (r"", t)dvdr'dv"dr"'dr"" 
+ J(VS® ■VS®)pW(r,t)dr + |(VS (1) - VS< 2 V 2) M)*} (31) 

In deriving these equations (Eq. p7|)-(|3TD) we have made use of the fact that 
for the ground-state = which implies that the current density p^ VS^ 

S(o_ 
dt 

use the first-order 

<9p« 



and the time-derivative ^757- vanish for the ground-state. In addition we also 



+ V • (p (0) VS (1) ) = (32) 



dt 

„<0 (t) = -^ + B g»(r,.) + \j J £^ J)P ^W (33) 
and the second-order 

^ + V • (p^VS^) + V ■ (p«VS«) = (34) 

" <2,(t) = + • vs(«) 

1 /■ ,«),., 



+ y MM)M?^p)r". t ) p(1 ' (r, - <)P '' ,(r "- tMr, ' ir " <35) 
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continuity and Euler equations obtained by expanding Eqs.(|2~I|) and 
Expressions for the average energies (Eqs.([2"7|)-(|3"T|)) clearly demonstrates the 
(2n + 1) rule of perturbation theory Energies up to order 3 are determined 
completely by p^ and S^. Similarly p^ and give energy up to the fifth- 
order. This is the (2n + 1) theorem of hydrodynamic perturbation theory 
in terms of the particle and the current densities. Moreover, even-order 
corollary of this theorem also holds true. Thus making E^ 2 ^ stationary with 
respect to and pW, leads to Eqs.(|32|) and d33|) . Similarly stationarity 
of E^ with respect to p^ and gives the correct perturbation equations 
(Eqs.(|34D and (|35D ) for p^ and S^ 2 \ The stationary nature of the even-order 



energies gives a variational method to obtain approximate solutions for the 
corresponding induced densities and currents. 

With this we complete the development of VP method in terms of particle 
and current densities in hydrodynamic formulation of TDDFT. In the next 
section we demonstrate the applicability of this theory by calculating the fre- 
quency dependent linear and nonlinear polarizabilities of inert gas atoms and 
comparing the results obtained with their wavefunctional counterparts. We 
then apply the formalism to calculate frequency dependent polarizabilities 
and plasmon frequencies of alkali metal clusters of large sizes. 

III. Application of hydro dynamical formalism 

To demonstrate the applicability of the formalism developed above we begin 
this section with the calculation of frequency dependent linear and nonlinear 
polarizabilities of inert gas atoms. In the present formulation we can calcu- 
late the nonlinear polarizabilities corresponding to the degenerate four wave 
mixing (DFWM) and DC-Kerr |35[ effect which are directly related to the 
fourth-order energy changes. On the other hand, unlike the orbital based 
Kohn-Sham approach, (2n + 1) theorem cannot be exploited to calculate 



p6|, [37], B3] the coefficients for the third-harmonic generation and electric 
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field induced second harmonic generation processes from only the second- 
order induced densities. In the following we will present the results for the 
nonlinear coefficients corresponding to the DFWM process only. 

As pointed out earlier, we perform our calculations using the variational 
property of the even-order energies. To this end we choose an appropriate 
variational form for the induced particle and current densities. For an applied 
potential of the form 

<lM) = <i(r)cos^ (36) 
with the spatial part given by 



v 



«(r) = £rcos£ (37) 



where £ is amplitude of the applied field, the time dependence of p and S at 
various orders can easily be inferred from Eqs(|32"l)-(l34T) as 



p (1) (r,t) = p (1) (r)coswt 

pW(r,t) = p( 2) (r)cos2^ + pJ 2) (r) (38) 

and 

S {1 \r,t) = S (1) {r)smut 

S (2) (r,i) = S (2) (r)sin2cut (39) 

Note that unlike the second-order particle density, the corresponding current 
has no constant term. This is consistent with the fact that the current arises 
due to the flow of electrons which causes the density to be time dependent. 
The spatial part of the induced particle and current densities are determined 
variationally by minimizing the appropriate even-order energies. For this 
purpose we choose the forms of the induced particle densities similar to the 



ones used previously |29|, ^0) for the calculation of static response properties. 
These are 

p (1) (r) = Ai(tW)cos0 (40) 
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and 



pf(r) = (A^(r) + A2(r)cos 2 ^)p (0) (r) + A 2 p(°)(r) 

p£\r) = (A°(r) + A°(r)cos 2 ^)p(°)(r) + Aop (0) (r) (41) 



with 



A,(r) = ^4^' ( 42 ) 
j 

where p'°'(r) is the ground-state density, a* are the variational parameters 
and As are fixed by the normalization condition for p( 2 \ To ensure satisfaction 

(2) (2) 

of the normalization condition at all times, p 2 and p are each normalized 
separately. These forms for the induced particle densities are motivated by 
the exact solutions [EM 011] for the hydrogen atom in a static field and have 



been shown [|29| , pOfl to lead to accurate static polarizabilities. On the basis 



of the continuity equation at each order and Eqs.(^0|)-(^), we choose 

S w (r) = A*(r)cos# (43) 

and 

S (2) (r) = (A*(r) + A{J(r) cos 2 0) , (44) 

where 

At(r)=J2by (45) 

3 

with b l j being the variational parameters to be determined by minimizing the 
average energy of respective orders. 

Application of hydrodynamical equations also requires approximating the 
KE and XC energy functionals. Based on our experience with the calculation 
of static linear and nonlinear polarizabilities we approximate them by their 



static forms. Thus for the KE, we use the von Weizsacker ETJ functional. 



TM = \j^y^ («) 
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For a discussion on the rationale behind choosing this functional, we refer the 
reader to the literature [^U], [3U], K|. For the exchange energy functional 
we use the adiabatic local-density approximation (ALDA) given by the Dirac 
exchange functional [33] 



EM 

a 



[v)dr 



(47) 



The correlation energy functional within the ALDA is represented by the 



Gunnarsson-Lundquist (GL) parametrization [45]. Thus 



EM = J e c (p)p(r)rfr 



with 



where 



e c (p) 



-0.0333 



1 + x 3 ) ln(l + -) + -x - x 2 
x 2 



x 



11.4 



(48) 



(49) 



(50) 



and r s = (J;-) 3 measures the radius in atomic units of a sphere which 
encloses one electron. 

The theory presented here treats the non-interacting KE exactly for sin- 
gle orbital systems (hydrogen and helium atoms). We have checked this by 
calculating the linear and nonlinear polarizabilities of H and He atoms. They 
match well with the corresponding wavefunctional results. The real test of 
the theory is therefore when it is applied to systems with more than one 
orbitals. We now present the results of these calculations by first discussing 
the frequency dependent polarizability numbers for the inert gas atoms of 
neon and argon. The ground-state electronic densities of these atoms are ob- 
tained by employing the van Leeuwen and Baerends (LB) ||6| potential. We 
use this potential as the orbitals generated by it have the correct asymptotic 



nature so that they lead to accurate values for the static [[47]] and frequency 
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dependent response properties. Here, instead of using the orbitals, we 



are using the ground-state densities generated by this potential. 

In Figs. 1 and 2 we show the linear polarizabilities for neon and 
argon atoms, respectively, as a function of lu, obtained by hydrodynamic 
approach. We compare these results (represented by open squares) with those 
obtained [|38| within the Kohn-Sham formalism shown in the figures by filled 
squares. As is evident, the frequency dependence matches quite well with 
the KS approach. This along with the zero frequency results demonstrate 
that the hydrodynamic theory is capable of giving reasonable estimate of 
dynamic polarizabilities in the optical range. Notice though that in the 
present approach the increase of a(u) with respect to u> is slightly less. 

To further quantify our results, we have fitted the frequency dependent 
polarizabilities with the formulae [[48| 



a(u) = a(0) (l + C 2 uj 2 ) (51) 

for small frequencies (up to u — 0.05 a.u.). In Table I we give the results for 
C*2 obtained from both the hydrodynamic and the orbital based calculations 



38| . As anticipated from the discussion above, the values of C 2 obtained 
from hydrodynamic formulation are close to but slightly smaller than their 
wavefunctional counterparts. 

Next, to study the performance of hydrodynamic approach in calculation 
of nonlinear response properties we calculate the coefficient corresponding 
to the DFWM phenomenon. These results are presented in Table II for 
two different frequencies, A = 10550A (a; « 0.0433 a.u.) and A = 6943A 
(u w 0.0657 a.u.). Here also we compare the present results with the cor- 
responding numbers obtained by the TD Kohn-Sham method. From 
this Table we again observe that the results for DFWM coefficients at both 
the frequencies are also quite close to, but lower than, the corresponding 
wavefunctional numbers. Notice that the numbers for hyperpolarizabilities 
are also underestimated by the hydrodynamic approach and the maximum 
deviation from the TD Kohn-Sham number is about 10%. However, with 
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the increase in frequency the difference between the hydro dynamical and the 
wavefunctional results will get larger. Nonetheless, the results obtained are 
reasonably accurate for frequencies up to uj ~ 0.05 a.u.. This is quite encour- 
aging since the experimental measurement of nonlinear coefficients fall within 
this frequency regime and also the numerical effort required for the hydrody- 
namical calculation is substantially less in comparison to the wavefunctional 
approach. 

Motivated by these results, in the next section we apply the hyrodynamic 
approach to calculate frequency dependent response properties and plasmon 
frequencies of metal clusters. As is well known, orbital based calculation for 
such systems are computationally demanding because of the large number of 
orbitals involved as the cluster size grows. 



Response properties of metal clusters 

The hydrodynamic approach developed above is particularly useful for sys- 
tems where an orbital based theory cannot be applied. Clusters are one such 
class of systems. These are made up of tens to thousands of atoms with 
properties distinct from the bulk properties of the constituent material. Fur- 
ther, various properties of clusters evolve in a well defined manner as their 



size grows. This was demonstrated by Knight et al. ||^, |5(|in their study of 
alkali metal clusters. Since then metal clusters have been studied [0, |53], |54fl 
quite extensively. One of the simplest model that describes average prop- 
erties of these systems correctly is the spherical jellium background model 
(SJBM) pH |53| , |54|. In small size clusters, Kohn-Sham LDA equations can 
be easily solved within this model. On the other hand, for large clusters one 
switches over to the density based theories the mostly applied one 

has been the extended Thomas- Fermi (ETF) theory |55| . In this paper also 
we use the density obtained by the ETF to calculate frequency dependent 
response properties (both linear and nonlinear) of alkali metal clusters with 
the number of atoms up to 5000. In the past, dynamic linear polarizabilities 
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of these clusters have been studied using the time-dependent Kohn-Sham 
theory pTI, 0, [5SJ but because of the difficulties mentioned above, the size 
up to which one could go has been limited. 

In the ETF method the ground-state density is obtained by minimizing 
the energy functional 

E[ P ] = T s ETF [p] + E H [p] + E L x ? A [p] + / Vj(v)p(v)dr + Ej (52) 

where V/ and Ej are the potential and the total electrostatic energy, respec- 
tively, of the ionic background. The functional Tf TF is the non-interacting 
KE functional included up to the fourth-order in density gradients. It is 
given as |59 



T s ETF [p]=T^[ P ]+TP[p]+T^[p] 



(53) 



where 

T (0) 

T (2) 

s 

T (4) 



(3n 2 ) 2 / 3 j pldr 
1 f|Vp| 2 



72 



-dr 



P 



540(3tt 2 



'vy 



9 (vy 

8 V P , 



, p , 



1 

3 



, p , 



Ejf[p] is the Hartree energy and for E^ is the LDA XC energy. For this we 



use the GL parametrization 



The energy above is minimized by taking 



the variational form [55] for the density to be 



P(r) = Po 



1 + exp 



R 



a 



(55) 



where R, a and 7 are the variational parameters and p$ is fixed by the 
normalization condition for each set of these parameters. The density so ob- 
tained gives results which are quite close [^] to the results of more accurate 
KS calculations of several properties. We use the ground-state densities ob- 
tained by this method as the input for the calculation of response properties. 
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We perform our calculations for sodium clusters with r s = 4.0, where r s is 
Wigner-Seitz radius of metal. Before presenting our results we point out that 
that the KE functional used to obtain the ground-state density and that for 
calculating the response properties are different. This is because whereas 
ETF functional is good for the total energies, it does not give the changes 
in the energies accurately. As mentioned earlier, for this purpose we use the 
von Weizsacker Jllj] functional. 

First we present the results for static linear polarizabilities. Although 
polarizability of metal clusters has been investigated extensively in the past 
PH, |57| , |58j, these studies have been restricted to clusters with number of 
atoms up to 200 because of the use of the orbitals in the calculations. We 
perform our study for clusters up to 5000 atoms. The variational forms for 
the induced particle and the current densities are chosen to be similar to 
those used in the atomic case. In Fig. 3 we show plot of static polarizability 
a(0) in the units of -Rjj as a function of Rq (where Rq = r s Nz, denotes the 
radius of cluster). It clearly shows that results of our calculation match quite 
well with the results obtained by Kohn-Sham approach for small (up to 196 
atoms) clusters. As the size of cluster grows the polarizability approaches 
the classical limit of Rq (that is — ► 1)- This is exhibited rather clearly 
in Fig. 3. 

Having obtained the static polarizabilities of alkali metal clusters accu- 
rately, we next study the dynamic response properties of metal clusters fo- 
cussing our attention particularly on the dipole resonance. The classical 
theory of dynamic polarizability predicts a single dipole resonance at the 
frequency given by |K| (in a.u.) 

^Mie = J —, (56) 
V ^ s 

which is equal to l/\/3~ times the bulk plasma frequency. The TDDFT results 
for the dipole resonance follow |53|, [54| the Mie result only in a qualitative 
way. The resonance peak corresponding to the photo absorption spectra of 
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these clusters exhaust about 70%-90% of the dipole sum rule and is red- 
shifted by about 10%-20% from the classical Mie formula 51]. 



In our work we estimate these resonance peaks from the frequency depen- 
dent polarizabilities by approximately locating the frequency at which a(u) 
becomes very large. These results are presented in Table III along with the 



results obtained by Brack |5j| using the RPA sum rules. It is clear from 
Table III that the dipole resonance frequencies of metal clusters obtained by 
hydrodynamical approach to TDDFT are quite accurate over the range of 
clusters studied. Further, the accuracy is better for larger clusters. We also 
find that with the increase in particle size the dipole resonance frequency 
approaches the classical Mie resonance frequency, which in the present case 
of r s = 4.0 is 0.125 a.u.. 

Next we discuss the results obtained for the nonlinear polarizabilities of 
metal clusters by the present approach. To the best of our knowledge hyper- 
polarizabilities for these systems have not been calculated before the present 
work. In these calculations we are restricted to clusters with maximum of 
only 300 atoms due to computational difficulties. Since electrons in metal- 
lic cluster are highly delocalized, we expect that the nonlinear response of 
these system should be quite large and increase rapidly with the size of the 
clusters. To ascertain how does the static hyperpolarizability 7(0) scale with 
the size of clusters, in Fig. 4 we plot 7(0) versus a(0) on a log-log scale. The 
line in this figure represents the best fit to hyperpolarizability versus polar- 
izability numbers obtained by us. It is seen from figure that for the clusters 
studied by us, the hyperpolarizability is linearly proportional to the linear 
polarizability. Since a(0) scales linearly with N, 7(0) also varies in the same 
manner. This is a surprising result since in the atomic case we have seen that 
7(0) increases much more rapidly than a(0) does. This could be because the 
electrons in metal clusters are much more mobile and therefore screen the 
applied potential very efficiently 

We have also studied the frequency dependent hyperpolarizabilities 7(0;; u, —u, u>) 
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and found it increasing with to. Variation of 7(0;; u, —uj, uo) with uo is shown 
in Fig. 5. It becomes quite large (by an order of magnitude in comparison 
to the static result) at approximately half the dipolar resonance frequency 
obtained from a(u) (Table III). This demonstrates the inherent consistency 
of the theory. 

Our study above has been done for clusters with number of atoms up 
to 300. However, the trends obtained should continue as the size grows. 
Slow increase of 7 with N is consistent with the fact that the classically 7 
for spherical metal particle is zero. One reason why computation becomes 
difficult for large clusters, we suspect, is that the variational forms chosen 
for the second-order particle and current densities may not be appropriate 
for very large clusters. As such investigations for larger clusters relegated to 
the future studies. 

Concluding remarks 

In this paper we have developed the time-dependent perturbation theory for 
periodic (in time) hamiltonian in terms of the particle and current densi- 
ties of electrons. For this we have employed the hydrodynamic equations of 
TDDFT. Application of the theory developed requires that the energy func- 
tional be approximated. We have demonstrated that with the von Weizsacker 
fil| functional for the KE and the ALDA for the XC energy, the theory leads 
to reasonably accurate results for dynamic response properties, both linear 
and nonlinear of atoms. Having established that we have applied the theory 
to study response properties of metallic clusters with number of atoms up to 
5000 within the SJBM. Of particular interest is how the hyperpolarizability 
varies with the size of these clusters. Although it is zero classically, our study 
shows that it increases linearly with the number of atoms in the cluster. 

Acknowledgement: We thank Prof. M. Brack for sending us his pro- 
gram to calculate ground-state density using the ETF approach. 
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Table Captions 

Table I C2 for inert gas atoms obtained by using hydrodynamical and wave- 
functional approaches. 

Table II DFWM coefficient 7(0;; 00, —ou, u) (in atomic units) for inert gas 
atoms using hydrodynamic approach. 

Table III Estimate of dipole resonance frequencies (in atomic units) of 
some metal clusters by using hydrodynamic approach. 
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Table I 



Atoms 


c 2 


c 2 




(hydrodynamical) 


(wavefunctional) w 


He 


1.12 


1.12 


Ne 


0.82 


1.04 


Ar 


2.16 


2.65 


Kr 


2.79 





(a) Ref. | 



Table II 



Atoms 


A = 


6943A 


A = 


10550A 


hydro dynamic 


wavefunctional^ 


hydrodynamic 


wavefunctional^ 


He 


44.47 


44.57 


43.50 


43.58 


Ne 


83.38 


94.06 


82.58 


91.65 


Ar 


1095.2 


1226.1 


1039.3 


1154.8 


Kr 


2392.7 




2229.5 




Xe 


5821.1 




5295.3 





(a) Ref. H 



Table III 



N 


Present 


RPA (aj 


8 


0.100 


0.113 


100 


0.105 


0.1198 


500 


0.1165 


0.1219 


1000 


0.121 


0.1226 


5000 


0.1223 


0.1236 



(a) Ref. [0 
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Figure Captions 

Fig.l Plot of a(uj)/a(0) as a function of frequency u for neon. The open 
and closed squares represent hydrodynamical and wavefunctional results 
respectively. 

Fig. 2 Plot of a(uj) / a(0) as a function of frequency u for argon. The open 
and closed squares represent hydrodynamical and wavefunctional results 
respectively. 

Fig. 3 Static polarizability a(0) in the units of R$ of alkali metal clusters 
as a function of Rq. The filled squares represents the results of Kohn-Sham 
calculations [p8"| . 

Fig. 4 Plot of 7(0) versus a(0) of alkali metal clusters. 

Fig. 5 Plot of 7(0;; u>, —uj, uj) in the units of R$ as a function of frequency 
uj for a metal cluster with 100 atoms. 
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